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1. INTRODUCTION 

Dynamic magnetic resonance imaging (DMRI) data is naturally stacked in two dimensions and 
involves entities that are evolving in time. Magnetic resonance imaging (MRI) is very important in many 
clinical examinations including vocal tract imaging, brain, cardiovascular, abdominal, and pulmonary [1], 
[2]. Any voluntary or involuntary motion generated by muscles or other physiological movements during the 
imaging procedure adds to the complexity of MRI image manipulation. Therefore, and because of the slow 
acquisition nature of DMRI, it usually compromises the quality of output images. In other words, for accurate 
quantification of myocardial perfusion in MRI data, it is essential to have an extended slice, a suitable spatio- 
temporal resolution, coverage, and a high signal to noise ratio (SNR) [3], [4]. 

To overcome aforementioned challenges, including parallel imaging, and their combination with k-t 
spatio-temporal models, many schemes were proposed in the literature [5]-[10]. However, these schemes are 
suitable for specific applications and tradeoff the quality with the complexity cost. Machine learning based 
schemes were proposed to solve the problem of motion estimation and deformation correction [11]-[14]. 
However, it does offer substantial evidence for the clinical significance of such a solution in the sense of 
motion estimation below the sampling rate. 

For some real-life applications, including medical imaging [15], the Nyquist rate is not applicable or 
viable, or even physically possible to build such a device, an explosive growth in the signal processing 
approach of compressed sensing (CS) for perfectly reconstructing data that are significantly sampled below 
the Nyquist rate [16]. Since the reconstruction of medical images required vectorization or matricization of 
the underlying data [17], CS is directly applied to multidimensional structure of dynamic MRI data. 
Dwivedi et al. [18] proposes a learning method using a tensor dictionary under an orthonormal constraint to 
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reconstruct DMRI as an alternates sparse coding algorithm. Lingala et al. [19] introduces a novel CS 
algorithm for jointly deformation correctness accelerated dynamic MRI, with a capability of handling a wide 
range of compactness/sparsity. Quantification feasibility for a whole-heart, non-electrocardiograph (ECG)- 
gated, cardiac-phase-resolved myocardial perfusion MRI framework and free-breathing was developed and 
tested in [20] using the CS framework. These methods have the ability to recover the data successfully for 
minimal inter frame motion, but they are sensitive to large inter frame motion. 
Lakshminarayana and Sarvagya [21] used CS in medical images with hybrid compression methods. 

In this piece of work, a novel estimation/correction approach were propose to reconstruct dynamic 
MRI images from under sampled measurements. The dynamic images and inter frame motion were jointly 
estimated. We used the variable splitting to decuple the main problem into three simpler sub-problems. Most 
importantly, we lowered the complexity cost of the deformation registration scheme (step), which registers 
each frame in the dataset to its correspondence step by step, through using patchMatch (PM) process instead. 
Finally, we illustrate the utility of our proposed framework in the context of myocardial perfusion. 

We organized the paper as follows: section 2 discusses the research method including the 
deformation-corrected CS framework and proposed algorithm. In section 3 we demonstrate the use of the 
proposed framework. Section 4 concludes our results. 


2. RESEARCH METHOD 

In this framework, it has been proposed a combining deformation-correction compressive sensing 
(DC-CS) approach using a global search based on PM novel random search algorithm. Under the 
presumption of exploiting redundant in similar patches across different subjects, it has been used NNF in 
association with the DC-CS approach for MRI image recovering. 


2.1. DC-CS for dynamic magnetic resonance imaging 

Inter frame motion recovery scheme can be modelled as low rank matrix problem which can be 
perfectly solved using compress sensing techniques. DC-CS minimizes the sensitivity of compressed data by 
estimating the dynamic images jointly with the inter frame motion. It was presumed to have a compact 
representation of the deformation corrected signal rather than the original signal. Similarly, we simplify the 
under-discussion problem using variable splitting into three smaller problems including a shrinkage-based 
denoising step, a quadratic optimization step, and deformable registration step. Also, in order to utilize the 
advantage of low rank structure in the construction, we considered various use of compactness priors 
including temporal finite difference domain, sparsity in the temporal Fourier domain, and nuclear norm 
penalty. 


2.2. Problem formulation 

The ultimate object of our framework is to recover the dynamic dataset f (x,t) from it’s under 
sampled Fourier noisy measures b(k;, tj) where x denotes the spatial variable and t refers to time. Also, we 
would like to simulate the result using complexity low-cost algorithm. The dynamic MRI measurement 
process can be modeled in vector space as follows (1): 


b=A(f)+n (1) 


where A is fourier transform evaluating operator from the sampling location in the set (kj, tj), j=1..s—1. 
Lingala et al. [19] proposes a deformation corrected dataset to overcome the limitation associated with 
dynamic MIR deformation correctness, Ig. It was assumed f to be compact-sparse and considerably 
structured. 

In this work we consider simultaneously recovery of the deformation parameters ®(x,t) and the 
dynamic images f (x, t) from under sampled data b(k, t) using the following general minimization scheme as (2): 


2 
0} = min AC) — dll, +48 Co. f) & 
f f such that mmo) 


Here A is the fourier sampling operator, Tọ is image operator for non-rigid warping; 0(f,t) is deformation 
parameters which represents pixel wise displacements due to motion. This parameter is estimated from the 
under sampled data. The proposed scheme also yields Tgcy+).f that is version of f, to represent the 
deformation corrected, as a by-product. @(u) is an arbitrary prior to exploit the redundancy in the data and 
finally © is the set of coefficients that represent the deformation field like B-spline. The authors used a 
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variable splitting approach to decouple the original problem to simpler sub-problems. An auxiliary variable g 
was introduced to split the prior from the deformation term which enables formulating the problem using the 
following constraints as (3): 


(fey = min HAC) — Pll, +480) & 
f,8.9 such that Tg, f =g9,0€ 90 

Figures 1(a)-(f) shows the DC algorithm, as a demonstration in the context of MRI myocardial 
perfusion. Hence, we can show some results on prospectively and retrospectively undersampled nuclear 
perfusion imaging (MPI) using radial trajectory. DC-CS reconstructions shows less motion artefacts and 
temporal blurring comparing to other CS schemes. Even though DC shows some promising results but still 
got some drawbacks. Some downsides of the algorithms are higher computational complexity, explicit 
motion compensation, slow running (mainly because of the deformation step, demon’s algorithm), high 
memory demand and didn’t use different choices of sampling patterns. To alleviate the drawback of the DC 
algorithm with confronted problem, we will use PM scheme. 


few image time frames Image time profile Temporal Fourier 


original 


deformation corrected 


(d) (e) (f) 


Figure 1. Free breathing multi-parametric MRI (MPMRI) with and without deformation correction 
(a) dynamic myocardial perfusion MRI frames, (b) the rippled time profile of the original datasets along the 
white dotted lines on the frames because of breathing, (c) temporal Fourier domain representation of the 
original dataset which shows the temporal harmonics corresponding to the motion are compensated in (f), (d) 
the deformated-corrected version of (a), (e) the free of ripples time profile of the deformated image along the 
white dotted lines on the frames, and (f) temporal Fourier domain representation of the deformation 
correction dataset which shows sparse representation 


The registration was, performed on the free breathing dataset itself. By using the demon’s 
registration algorithm, starting from the second frame, the deformations were found by matching the n-th 
frame in the moving, sequence to the (n+1)-th of the deformed scene. Solving the optimization problem by 
relaxing the unconstrained problem yields (4): 


2 b 2 
(F,0)° = min A — Al, +A((g) +>|Ile,f - gll) A 
f.g such that 0 E€ © 
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Beta is the penalty parameter enforces the constraint. Using the iterative minimization scheme to 
solve the problem above by solving one variable and fixing the rest. They reiterate the three sub-problems 
until they converge. 

- g sub-problem (spatio-temporal denoising/dealiasing to minimize residual motion as (5)-(7): 


g` = min O(g) + IIo, f — gll where q = To, f (5) 
a = 2a a a 
g = min=|9lli, + llâ- 93 (6) 
g 
— f(a) 2 
g = EL (IF @I-=) 7) 


- f sub-problem (reconstruction update) as (8): 

min||A(f) ~ Bll, + Ifo. f = 911 (8) 
- Theta sub-problem (motion estimation) as (9): 

min |ITo,f — gll (9) 


The dynamic f (x,t) is registered frame by frame with a reference scene g(x,t). Here, an iterative 
demon’s algorithm was used to approximate the registration subproblem. As demon’s is an iterative 
algorithm, the displacement field parameter 6 (x, t; ) is updated in each iteration j by adding a force field. 

Factor Un, On+1 = On + Un, Where u, evaluated in the previous iteration (n — th iteration). The force 
field is derived from intensity differences in the reference and target images. The risk of converging into a 
local minimum is still exist, therefore; the author employs continuation strategies to minimize its risk 
initializing beta parameter with very low values and increase it gradually. To speed up the convergence of the 
algorithm and to reduce local minima effects, the demon’s implementation recommends a continuation on a 
force strength parameter. The resulting algorithm decreases the presence of local minima convergence 
through employing efficient continuation strategies. DC-CS considers a deformable registration scheme that 
registers the frames from the existing dataset to their corresponding. 


2.3. PatchMatch technique for structural image editing 

The work in [22] considers one of the landmarks for all time in image editing for its interactive 
ability of photo manipulation using fast randomize algorithm. The key strengthen of the framework is its new 
nearest-neighbour field (NNF) search algorithm which uses random sampling to find some good patch 
matches which quickly propagate to the surrounding areas through the cohesive nature of the imagery. The 
previous NNF algorithms attain high synthesis quality of patch optimization, but with a trade-off with the 
complexity which the bottleneck of such algorithm. In contrast, the technique in [22] not only achieves high 
quality results, but also provides almost real time photo editing tools for image application including image 
retargeting, completion, and reshuffling. 

Most of the other algorithms which try to accomplish image reconstruction look in the image for 
regions that are like the one that was removed and borrow some information to fill out the process. Unlike these 
algorithms, when we are trying to do another lookup process, our algorithm does not restart the patch matching 
process. However, when good correspondence is found, it just searches nearby because that is where we can 
find the most useful information. The main idea here is that exploiting the local coherence in the synthesis 
process, speeds up the filling process through minimizing the search space to the neighbours of the patch local 
location. We borrow the coherence propagation assumption of our algorithm from the work in the [23]. 

In the context of the image editing application, our NNF algorithm, which is illustrated in Figure 2, 
has three phases: i) initially, patches possess random assignments; ii) the blue patch checks above/green and 
left/red neighbours to check whether there is a possible improvement to the blue mapping, and propagating a 
possible good match; and iii) the patch randomly searches the concentric neighbourhoods for improvements. 
In other words, a random offset (sometimes prior information) will be filling the NNF to apply an iterative 
update process. Adjacent pixels are occupied by good patch offsets. PM can determine an approximate 
‘nearest neighbour’ match that could conceivably replace a given patch in the image. Once a reasonable match 
is found, it tries to improve this patch by randomly searching for better matches around the target position. 
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(a) Initialization (b) Propagation (c) Search 


Figure 2. Scheme phases using a randomized approximation algorithm [23] 


2.4. Proposed patchMatch for DC-CS scheme magnetic resonance imaging recovery 

One of the main contributions of this work is introducing a global search strategy instead of the DC 
registration step, guaranteeing a non-explicit motion estimation that avoids any spatial alignment or 
registration of the images; in addition to use PM. Scheme to decrease the computational cost to the minimum. 
We propose a framework to combine the DC-CS approach using a global search based on the PM algorithm. 
We assume that information redundancy in similar patches across different subjects can be exploited. The 
basic patch match algorithm allows for the best matches between a source and target image to be found 
quickly. Most importantly we can change the slow and complex registration step by using the patch-based 
algorithm to implicitly estimate motion and do image reconstruction. 


3. RESULTS AND DISCUSSION 

The main reason behind thinking that our proposed scheme can be utilized in the image recon is that 
this patch-based method uses NNF which is good for motion estimation and the randomized search followed 
by propagation is not just fast, but easier to implement smoothness enforcement so we make sure the 
computational efficiency is kept lower. The basic patch match algorithm allows for the best matches between 
a source and target image to be found quickly. Most importantly we can change the slow and complex 
registration step by step using the patch-based algorithm to implicitly estimate motion and do image 
reconstruction. We did the simulation in Figure 3 with and without using the propagation step. We found that 
the result using the propagation step is pretty smooth and converges even better. 


x 
Propogation On 
— Propagation Off 


Figure 3. Error with (smooth) and without (zigzag) using propagation step 


Figures 4 shows the implementation of PM which is done in C** with a MEX interface to 
MATLAB. The starting points for them were chosen to be zeros and the neighborhood size was set of 20. 5 
patches searches are involved and with a size of 5*5 each, thus we may recover five images from the original 
image. We have done the experiments using 256*256 vase image and 190*90 MPI. One observation is that 
PM preserves the edges well and gives good results. The 1“ image in both figures always looks the best as it 
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is simply the translational version of the patches of the original image while the rest four images are obtained 
from the four neighboring patches within the search vicinity. 

As PM breaks down an image into a series of small, overlapping rectangles of fixed size; it replaces 
the wrapping operator in DC by this scheme. More specifically, we replace the forward operator by 2" 
neighbourhood PM and the backward one as the 1* neighbourhood which is simply the translated version of 
the image. We can also save time by doing PM only once by saving the image as an array of indices before 
applying the PM so we retrieve it easily by applying the inverse wrapping operator. By doing so we can make 
sure that the patches capture the motion in the frames without any complicated or explicit motion estimation. 

Figure 5 shows one frame of the original image and the recon using the PM recovered from an 
acceleration factor of 5 using the radial undersampling operator. The shrinkage and f calculation steps were 
kept as they are in DC. We can see that the recovered image is almost an artefact-free and crispy. Minimizing 
the input and output number of iterations set in the DC algorithm might speed up the whole process. It is 
worth to mention that there are many other algorithms that can be tried to use which might improve the 
quality and the complexity of the work. We can also collaborate this algorithm with a neural network which 
can be a future improvement over this work [24]-[27]. 


original SNR=15.2 


Figure 5. One frame of the original image and the recon using the PM recovered from acceleration factor of 5 
using radial undersmapling operator 
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4. CONCLUSION 

In this framework, a novel estimation/correction approach was proposed for dynamic MRI image 
reconstruction from its under sampled measurements. The dynamic images and inter frame motion were 
jointly estimated. Also, we proposed a PM algorithm to be used instead of DC algorithms which has a high 
complexity that is the main challenge of data reconstruction. We used a global search strategy instead of the 
DC registration step to guaranteeing a non-explicit motion estimation that avoids any spatial alignment or 
registration of the images. One of the limitations we have faced is that the PM algorithm accepts only real- 
valued images that we are planning to address in the future. 
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